use "Samuels_Thomson_Brazil_Tractors_Land_Invasions.dta", clear
set scheme plotplain

*********************
* Describe the data *
*********************
*** The Invasions Data
hist albertus_88_90 if year==1985, ///
discrete percent xsize(4) ysize(4) ///
title("1988-90") xtitle("Invasion Count")
graph save invasion_hist_88_90, replace

hist albertus_96_98 if year==1996, ///
discrete percent xsize(4) ysize(4) ///
title("1996-98") xtitle("Invasion Count")
graph save invasion_hist_96_98, replace

hist albertus_06_08 if year==2006, ///
discrete percent xsize(4) ysize(4) ///
title("2006-08") xtitle("Invasion Count")
graph save invasion_hist_06_08, replace

hist occs_17_19 if year==2017, ///
discrete percent xsize(4) ysize(4) ///
title("2017-19") xtitle("Invasion Count")
graph save invasion_hist_17_19, replace

graph combine invasion_hist_88_90.gph invasion_hist_96_98.gph ///
invasion_hist_06_08.gph invasion_hist_17_19.gph, ycommon cols(2) ysize(4) xsize(4)


*** Correlation tables

corr binary_88_90 std_tractors std_income std_landgini std_rural if year==1985

corr binary_96_98 std_tractors tri_mean_abs std_income std_landgini std_rural if year==1996

corr binary_06_08 std_tractors tri_mean_abs std_income std_landgini std_rural if year==2006

corr binary_17_19 std_tractors tri_mean_abs std_income std_landgini std_rural if year==2017

*** The Tractors data
hist tractors if year==1985, ///
discrete percent xsize(4) ysize(4) ///
title("1985") xtitle("Tractors")
graph save tractors_hist_85, replace

hist tractors if year==1996, ///
discrete percent xsize(4) ysize(4) ///
title("1996") xtitle("Tractors")
graph save tractors_hist_96, replace

hist tractors if year==2006, ///
discrete percent xsize(4) ysize(4) ///
title("2006") xtitle("Tractors")
graph save tractors_hist_06, replace

hist tractors if year==2017, ///
discrete percent xsize(4) ysize(4) ///
title("2017") xtitle("Tractors")
graph save tractors_hist_17, replace

graph combine tractors_hist_85.gph tractors_hist_96.gph ///
tractors_hist_06.gph tractors_hist_17.gph, cols(2) ysize(4) xsize(4)






************
* ANALYSIS *
************
sort codeibge year


*********
* Logit *
*********

*** 1988-1990

* Baseline, Table 1
logit binary_88_90 std_tractors i.state, vce(cluster state)
est store logit_88_90_baseline

*Main model, Table 1
logit binary_88_90 std_tractors std_landgini std_rural std_income std_area std_distance i.state, vce(cluster state)
est store logit_88_90_1

* Margins for tractors
margins, dydx(std_tractors) at(std_tractors = (-1(0.25)1)) predict(pr)
marginsplot, title("1988-90") ///
ytitle("Invasion Probability") ///
ciopts(recast(rarea) color(%10)) ///
addplot(kdensity std_tractors if e(sample) & std_tractors>=-1 & std_tractors <=1, yaxis(2) yscale(r(0 5) axis(2) off) lpattern(longdash)) ///
xsize(5) ysize(5) xlabel(#5) legend(off) ///
saving(logit_88_90.gph, replace)



*** 1996-1998

* Baseline, Table 1
logit binary_96_98 std_tractors invasions_lag i.state, vce(cluster state)
est store logit_96_98_baseline

* Main model, Table 1
logit binary_96_98 std_tractors std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_96_98_1

* Margins for tractors
margins, dydx(std_tractors) at(std_tractors = (-1(0.25)1)) predict(pr)
marginsplot, title("1996-98") ///
ytitle("Invasion Probability") ///
ciopts(recast(rarea) color(%10)) ///
addplot(kdensity std_tractors if e(sample) & std_tractors>=-1 & std_tractors <=1, yaxis(2) yscale(r(0 3) axis(2) off) lpattern(longdash)) ///
xsize(5) ysize(5) xlabel(#5) legend(off) ///
saving(logit_96_98.gph, replace)

* Margins for Land Gini
margins, dydx(std_landgini) at(std_landgini = (-1(0.25)1)) predict(pr)

***** Table 2
* Pesticide/fertilizer
logit binary_96_98 std_tractors std_pesticide std_fertilizer std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_96_98_2
* Albertus: Nearby Reform
logit binary_96_98 std_tractors ref_100km std_pesticide std_fertilizer std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_96_98_3
* Albertus: Include the interaction (do not report)
logit binary_96_98 std_tractors c.ref_100km##c.std_landgini std_pesticide std_fertilizer std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
** Hidalgo et al.: Ag Income
logit binary_96_98 std_tractors std_agincome ref_100km std_pesticide std_fertilizer std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_96_98_4
** Hidalgo et al.: Rainfall
logit binary_96_98 std_tractors rainfall_96_98 ref_100km std_pesticide std_fertilizer std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_96_98_5
*** Falcone Rosenberg: Crop Yields
logit binary_96_98 std_tractors rainfall_96_98 ref_100km std_pesticide std_fertilizer std_soy std_maize std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_96_98_6
*** F&R Interact with Tractors (do not report)
logit binary_96_98 c.std_tractors##c.std_soy rainfall_96_98 ref_100km std_pesticide std_fertilizer std_maize std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
logit binary_96_98 c.std_tractors##c.std_maize rainfall_96_98 ref_100km std_pesticide std_fertilizer std_soy std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)


** Robustness Tests
* 5-year window
logit binary_96_00 std_tractors std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_96_00_1
* 10-year window
logit binary_96_05 std_tractors std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_96_05_1



*** 2006-2008

* Baseline
logit binary_06_08 std_tractors invasions_lag i.state, vce(cluster state)
est store logit_06_08_baseline

* Main model, Table 1
logit binary_06_08 std_tractors std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_06_08_1

* Margins for tractors
margins, dydx(std_tractors) at(std_tractors = (-1(0.25)1)) predict(pr)
marginsplot, title("2006-08") ///
ytitle("Invasion Probability") ///
ciopts(recast(rarea) color(%10)) ///
addplot(kdensity std_tractors if e(sample) & std_tractors>=-1 & std_tractors <=1, yaxis(2) yscale(r(0 5) axis(2) off) lpattern(longdash)) ///
xsize(5) ysize(5) xlabel(#5) legend(off) ///
saving(logit_06_08.gph, replace)

* Margins for land Gini
margins, dydx(std_landgini) at(std_landgini = (-1(0.25)1)) predict(pr)

* Table 2
* Pesticide/fertilizer
logit binary_06_08 std_tractors std_fertilizer std_pesticide std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_06_08_2
* Albertus nearby reform variable
logit binary_06_08 std_tractors ref_100km std_pesticide std_fertilizer std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_06_08_3
* Albertus: Interact with land Gini (do not report)
logit binary_06_08 std_tractors c.ref_100km##c.std_landgini std_pesticide std_fertilizer std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
** Hidalgo et al.: Ag Income
logit binary_06_08 std_tractors std_agincome ref_100km std_pesticide std_fertilizer std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_06_08_4
** Hidalgo et al.: Rainfall
logit binary_06_08 std_tractors rainfall_06_08 ref_100km std_pesticide std_fertilizer std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_06_08_5
*** Falcone Rosenberg: Crop Yields
logit binary_06_08 std_tractors rainfall_06_08 ref_100km std_pesticide std_fertilizer std_soy std_maize std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_06_08_6


** Robustness Tests
* 5-year window
logit binary_06_10 std_tractors std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_06_10_1
* 10-year window
logit binary_06_15 std_tractors std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_06_15_1



*** 2017-2019

* Baseline
logit binary_17_19 std_tractors invasions_lag i.state, vce(cluster state)
est store logit_17_19_baseline

* Main model, Table 1
logit binary_17_19 std_tractors std_landgini std_rural std_income std_area std_distance invasions_lag i.state, vce(cluster state)
est store logit_17_19_1

margins, dydx(std_tractors) at(std_tractors = (-1(0.25)1)) predict(pr)
marginsplot, title("2017-19") ///
ytitle("Invasion Probability") ///
ciopts(recast(rarea) color(%10)) ///
addplot(kdensity std_tractors if e(sample) & std_tractors>=-1 & std_tractors <=1, yaxis(2) yscale(r(0 5) axis(2) off) lpattern(longdash)) ///
xsize(5) ysize(5) xlabel(#5) legend(off) ///
saving(logit_17_19.gph, replace)

*** Results Graph for naive Logit
graph combine logit_88_90.gph logit_96_98.gph logit_06_08.gph logit_17_19.gph, ///
rows(2) ycommon title() xsize(5) ysize(5)




*** MAIN RESULTS TABLES ***


***********
* Table 1 *
***********
esttab logit_88_90_baseline logit_96_98_baseline logit_06_08_baseline logit_17_19_baseline ///
logit_88_90_1 logit_96_98_1 logit_06_08_1 logit_17_19_1, ///
keep(std_tractors std_landgini std_rural std_income std_area std_distance) ///
order(std_tractors std_landgini std_rural std_income std_area std_distance)  ///
mtitles("88-90" "96-98" "06-08" "17-19" "88-90" "96-98" "06-08" "17-19") ///
b(2) aic se compress nogaps label obslast

***********
* Table 2 *
***********
esttab logit_96_98_2 logit_06_08_2  logit_96_98_3 logit_06_08_3 logit_96_98_4 logit_06_08_4 logit_96_98_5 logit_06_08_5 logit_96_98_6 logit_06_08_6, ///
keep(std_tractors std_fertilizer std_pesticide ref_100km std_agincome rainfall_96_98 rainfall_06_08 std_soy std_maize std_landgini std_rural std_income) ///
order(std_tractors std_fertilizer std_pesticide ref_100km std_agincome rainfall_96_98 rainfall_06_08 std_soy std_maize std_landgini std_rural std_income)  ///
mtitles("96-98" "06-08" "96-98" "06-08" "96-98" "06-08" "96-98" "06-08" "96-98" "06-08") ///
b(2) aic se compress nogaps label obslast ///
addnotes("Standard errors clustered by state." ///
"All models control for lagged invasions, municipality area and distance from capital.")




*** APPENDIX / ROBUSTNESS TESTS ***

***********************
* IV Analysis: Probit *
***********************

*** 1988-1990
* Main model, Table 1
ivprobit binary_88_90 (std_tractors = tri_mean_abs) std_landgini std_rural std_income std_area std_distance i.state, vce(robust)
est store ivprobit_88_90_1

margins, dydx(std_tractors) at(std_tractors = (-1(0.25)1)) predict(pr)
marginsplot, title("1988-90") ///
ytitle("Invasion Probability") ///
ciopts(recast(rarea) color(%10)) ///
addplot(kdensity std_tractors if e(sample) & std_tractors>=-1 & std_tractors<=1, yaxis(2) yscale(r(0 5) axis(2) off) lpattern(longdash)) ///
xsize(5) ysize(5) xlabel(#5) legend(off) ///
saving(ivprobit_88_90.gph, replace)


*** 1996-1998
* Main model, Table 1
ivprobit binary_96_98 (std_tractors = tri_mean_abs) std_landgini std_rural std_income std_area std_distance i.state, vce(robust)
est store ivprobit_96_98_1

* Margins for tractors
margins, dydx(std_tractors) at(std_tractors = (-1(0.25)1)) predict(pr)
marginsplot, title("1996-98") ///
ytitle("Invasion Probability") ///
ciopts(recast(rarea) color(%10)) ///
addplot(kdensity std_tractors if e(sample) & std_tractors>=-1 & std_tractors<=1, yaxis(2) yscale(r(0 5) axis(2) off) lpattern(longdash)) ///
xsize(5) ysize(5) xlabel(#5) legend(off) ///
saving(ivprobit_96_98.gph, replace)


*** 2006-2008
ivprobit binary_06_08 (std_tractors = tri_mean_abs) std_landgini std_rural std_income std_area std_distance i.state, vce(robust)
est store ivprobit_06_08_1

* Margins for tractors
margins, dydx(std_tractors) at(std_tractors = (-1(0.25)1)) predict(pr)
marginsplot, title("2006-08") ///
ytitle("Invasion Probability") ///
ciopts(recast(rarea) color(%10)) ///
addplot(kdensity std_tractors if e(sample) & std_tractors>=-1 & std_tractors<=1, yaxis(2) yscale(r(0 5) axis(2) off) lpattern(longdash)) ///
xsize(5) ysize(5) xlabel(#5) legend(off) ///
saving(ivprobit_06_08.gph, replace)


* 2017-2019
ivprobit binary_17_19 (std_tractors = tri_mean_abs) std_landgini std_rural std_income std_area std_distance i.state, vce(robust)
est store ivprobit_17_19_1

margins, dydx(std_tractors) at(std_tractors = (-1(0.25)1)) predict(pr)
marginsplot, title("2017-19") ///
ytitle("Invasion Probability") ///
ciopts(recast(rarea) color(%10)) ///
addplot(kdensity std_tractors if e(sample) & std_tractors>=-1 & std_tractors <=1, yaxis(2) yscale(r(0 5) axis(2) off) lpattern(longdash)) ///
xsize(5) ysize(5) xlabel(#5) legend(off) ///
saving(ivprobit_17_19.gph, replace)


***************************
* Tests of the Instrument *
***************************
* Simple bivariate regression
reg std_tractors tri_mean_abs if year==1985
est store ivtest_85_1

reg std_tractors tri_mean_abs if year==1996
est store ivtest_96_1

reg std_tractors tri_mean_abs if year==2006
est store ivtest_06_1

reg std_tractors tri_mean_abs if year==2017
est store ivtest_17_1

esttab ivtest_85_1 ivtest_96_1 ivtest_06_1 ivtest_17_1, ///
mtitles("1985" "1996" "2006" "2017") ///
b(2) se compress nogaps label obslast scalars(F)

* Look at relationship with fertilizer, pesticide
corr tri_mean_abs std_tractors std_fertilizer std_pesticide if year==1996

corr tri_mean_abs std_tractors std_fertilizer std_pesticide if year==2006

**************************
* Results of first stage *
**************************
esttab ivprobit_88_90_1 ivprobit_96_98_1 ivprobit_06_08_1 ivprobit_17_19_1, ///
keep(std_tractors:tri_mean_abs std_tractors:std_landgini std_tractors:std_rural std_tractors:std_income std_tractors:std_area std_tractors:std_distance) ///
order(std_tractors:tri_mean_abs std_tractors:std_landgini std_tractors:std_rural std_tractors:std_income std_tractors:std_area std_tractors:std_distance) ///
mtitles("1988-90" "1996-98" "2006-08" "2017-19") ///
b(2) se compress nogaps label obslast


**********************************
* IV 2nd stage for Appendix only *
**********************************
esttab ivprobit_88_90_1 ivprobit_96_98_1 ivprobit_06_08_1 ivprobit_17_19_1, ///
keep(std_tractors std_landgini std_rural std_income std_area std_distance) ///
order(std_tractors std_landgini std_rural std_income std_area std_distance)  ///
mtitles("88-90" "96-98" "06-08" "17-19") ///
b(2) aic se compress nogaps label obslast


** Graph of IVProbit Results
graph combine ivprobit_88_90.gph ivprobit_96_98.gph ivprobit_06_08.gph ivprobit_17_19.gph, ///
rows(2) ycommon title() xsize(5) ysize(5)


**************************************
* Robustness tests for Appendix only *
**************************************
esttab logit_96_00_1 logit_96_05_1 logit_06_10_1 logit_06_15_1, ///
keep(std_tractors invasions_lag std_landgini std_rural std_income std_area std_distance) ///
order(std_tractors invasions_lag std_landgini std_rural std_income std_area std_distance)  ///
mtitles("96-00" "96-05" "06-10" "06-15") ///
b(2) aic se compress nogaps label obslast ///
addnotes("Standard errors clustered by state." ///
"All models control for lagged invasions, municipality area and distance from capital.")